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— ■ ABSTRACT 

O' 

o: 

CN ■ We suggest a new approach that could be used for modeling both the large 

^ ! scale behavior of astrophysical jets and the magnetically dominated explosions in 

^ ! astrophysics. We describe a method for modeling the injection of magnetic fields 

I ! and their subsequent evolution in a regime where the free energy is magnetically 

CN ! dominated. The injected magnetic fields, along with their associated currents, 

,_H '. have both poloidal and toroidal components, and they are not force free. The dy- 

^ ! namic expansion driven by the Lorentz force of the injected fields is studied using 

\0 '. 3-dimensional ideal magnetohydrodynamic simulations. The generic behavior of 

^ ! magnetic field expansion, the interactions with the background medium, and the 

O ! dependence on various parameters are investigated. 

o| 

r-| ! Subject headings: magnetic fields — galaxies: active — galaxies: jets — methods: 

^', numerical — MHD 

^: 

^. : 1. INTRODUCTION 

>. 

r> I It is generally thought that energies (including magnetic fields) in jets and lobes of radio 

a I galaxies observed on large scales (from kpc to Mpc) are supplied from the central black hole 

accretion region, despite our poor understanding on how exactly this has occurred. Near the 
black hole, studies have focused on the initial formation of jets, in terms of jet coUimation 
and acceleration. Here, both hydromagnetic limit (e.g., Blandford & Payne 1982; Guyed & 
Pudritz 1997; Ustyugova et al. 1999; and others) and Poynting flux limit (e.g., Blandford 
1976; Lovelace 1976; Lynden-Bell 1996; Ustyugova et al. 2000; Li et al. 2001; Lovelace et al. 
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2002; and others) have been investigated, with various degree of success [see Ferrari (1998) 
for a relatively recent review]. 

On much larger scales, jet propagation, morphologies (such as lobe formation, hot-spots, 
filaments, etc.) and associated instabilities have been investigated extensively. In particular, 
much progress has been made in the so-called kinetic energy dominated regime, for both 
non-relativistic and relativistic flows. The pioneering work by Norman et al. (1988) set the 
stage for nearly two decades of studies on how highly supersonic, light (jet material density 
less than the surrounding density), and sometimes magnetized, collimated flows propagate 
in various (uniform and/or stratified) background medium (e.g., Norman et al. 1983; Marti 
et al. 1994; Clarke 1996; Bodo et al. 1998; Tregillis et al. 2004; and many others). Even 
though magnetic fields have been actively included in some of these studies, the overall 
behavior has been predominantly determined in the hydrodynamic limit (but see a recent 
paper by Nakamura & Meier 2004 in both hydro and Poynting dominated limits). 

The relative importance of kinetic energy and magnetic energy in jets and radio lobes 
has been a subject for extensive studies and debates. Better radio observations on large, 
> tens of kpc, scales (e.g., Owen et al. 2000) and recent X-ray observations of radio lobes 
of FanarofF- Riley type II (FR II) sources (e.g., Croston et al. 2005) have suggested that 
magnetic energies inside the radio lobes are significant (see also Kronberg et al. 2001). These 
observations bring out the questions whether jets/lobes could be magnetically dominated and 
what are the physical processes govern such flows on these large scales. The studies on this 
Poynting flux dominated limit in large scales have been relatively sparse. 

The purpose of the present paper is to present an approach to model the large scale 
magnetic structure generated by a purely magnetic energy input. We describe a method 
of injecting magnetic fields into a small volume and following the subsequent evolution of 
the injected magnetic free energy. This free energy causes magnetic fields to expand under 
their own Lorentz force. This expansion is subsequently affected by interacting with the 
background environments. We believe that this approach could in principle be applied to 
different astrophysical systems. Here, we concentrate on its applications to modeling the 
large scale morphologies of radio jets and lobes. In §2, we describe our basic assumptions 
and approach on how to model magnetic field injection. Numerical methods used in our 
simulations are also discussed. In §3, we present the simulation results on the evolution of 
magnetic fields for different parameters and setups. In §4, we summarize our results. 
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2. Basic Approach and Model Assumptions 

Since it is currently impossible to simulate the whole process from the jet formation near 
the black hole accretion disk out to radio lobes because of the vast length scale separation, we 
adopt the following general framework which separates the jet/lobe structure into different 
regions: I: the engine region, where the black hole accretion disk system produces magnetic 
fields and these fields are ejected above and below the disk. II: the propagation region, where 
energy, mass, and magnetic fields flow through the background medium. Ill: the termination 
region, where these flows are significantly slowed down or come to a stop. In our approach, 
for region I, we emphasize that we are not modeling the black hole accretion disk system and 
how magnetic fields are ejected from such systems. Instead, we attempt to model the likely 
outcome of the energy/mass supply from the black hole accretion system as having particular 
functional forms in non-force-free magnetic fields, along with certain injection rates. In this 
sense, we are focusing on the energy/mass flows in an external medium, and how they would 
eventually terminate (i.e. Regions II and III). In particular, we study the limit that the 
injected free energy is predominantly in magnetic fields, which is different from the previous 
approach where the flow is dominated by the kinetic energy (i.e., the injected kinetic energy 
is much larger than the injected magnetic and thermal energy). In the following sections, 
we will discuss the functional forms of injected magnetic fields and the numerical setups for 
modeling the nonlinear evolution of these magnetic fields. 



2.1. Basic Equations, Magnetic Fields, Mass Injections and Numerical Codes 

We solve the ideal MHD equations in three-dimensional Cartesian coordinate system 
{x,y,z]: 



dp 

d{pv) 



^^+V-(pv) = pinj (1) 



^^ , V ■ (pvv + Pg + Pb - BB) = (2) 

dE 

_ + V.[(E + Pg + PB)v-B(v.B)] = Ei,j (3) 

— -Vx(vxB) = Bi,j , (4) 

where all variables have their usual meaning. We use an ideal gas equation of state with an 
adiabatic index 7 = 5/3, which relates the internal energy density e and the gas pressure 
Pg as e = Pg/(7 — !)• The magnetic energy density is expressed as Pb = P^/2. The total 
energy E = pv'^/2 + Pg + Pb. We have neglected external forces such as gravity in these 
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simulations. This is not correct for some astrophysical environments such as jets in galaxy 
clusters. The main focus of this paper is to present the salient features of this new approach, 
and we will present more detailed astrophysical applications in future publications. 

We have added an injection term Binj in the induction equation to incorporate the 
external injection of magnetic fields (V ■ Binj = 0). The injection process is modeled as: 

• _ / 76 Bi„j for t < tinj , . 

^"'^ ~\ for t > ti,j ^^' 

where Binj describes the magnetic fields whose exact form is given in the following section. 
The 7fe(t) specifies the injection rate for magnetic fields, and it is usually taken as a constant 
over a finite time internal t < ti^j and as zero (e.g., the injection is turned off) for t > tinj. 
More general forms of injection can be easily included as well. 

In addition, we have included mass injection in the central region, motivated by the pos- 
sibility that matter can be ejected together with magnetic fields away from the central source. 
Note that both the magnetic field and mass injections violate the ideal MHD condition, in 
the sense that the total mass on magnetic field lines are allowed to change. Furthermore, 
adding mass unto field lines helps to alleviate the numerical problems of having extreme low 
density regions after magnetic fields have undergone large expansions. Without the detailed 
knowledge on how and what amount of mass could get unto the field lines, we adopt the 
following simple formula: 

Anj = IpPo exp[-{r^ + z'^)/rl] , (6) 

where Vp and 7^ are the characteristic radius and rate for mass injection. The injected mass 
is further assumed to have the same velocity and temperature as that of the gas in the 
injection region. 

Furthermore, there is an energy injection along with the magnetic field and mass injec- 
tions. The energy injection rate from the magnetic fields injection can be expressed as 

-^inj = B ■ Binj = 7bBinj ■ B , (7) 

while in general both B and Binj have temporal and spatial dependences. Note that there is 
a small amount of thermal and kinetic energy injected, accompanied with the mass injection. 
But this energy is negligibly small compared with the injected magnetic energy. 

All simulations are performed using a new ideal MHD package developed at Los Alamos 
(Li & Li 2003). This package uses high-order Godunov-type finite volume numerical methods. 
These methods conservatively update the zone-averaged fluid and magnetic field quantities 
based on estimated advective fluxes of mass, momentum, energy and magnetic field at zone 
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interfaces. Due to the discontinuities (contacts and shocks) in the solutions, Roe's Riemann 
solver with an entropy-fix (Powell et al. 1994) is used to calculate these fluxes. It usually 
solves the total energy equation as the ideal MHD equations are written into the conservative 
form. To maintain pressure positivity in regions where the pressure can become very small 
or negative, we also solve two auxiliary equations: the internal energy equation and the 
modifled entropy equation, 

_ + V-(ev) = -PgV-v (8) 

— + V-(5v) = (9) 

where S = Pg/p*^^~^^ is deflned as the modifled entropy. The exact usage of these two 
equations is described in Ryu et al. (1993) and Balsara & Spicer (1999a). 

The divergence free condition of the magnetic fleld V ■ B = is ensured by a constrained 
transport (CT) scheme flux-CT (Balsara & Spicer 1999b). The CT scheme can yield mag- 
netic flelds that are different from the updated values based on the Godunov method. The 
difference may be large enough to lead to negative pressure. To minimize the change in en- 
ergy conservation, we update the total energy with the newly-calculated magnetic flelds only 
when the pressure becomes negative. The whole package is parallelized via message-passing 
interface (MPI). A typical run made on parallel Linux clusters with 64 processors takes 3 
hrs. 



2.2. Injected Magnetic Field Configuration 

We determine the injected magnetic flelds using three key quantities: the length scale 
of the injection region (which we designate as r = 1), the amount of poloidal flux (\E'p), and 
the poloidal current (J^). For simplicity, we also assume that the injected magnetic flelds 
Binj are axisymmetric. In cylindrical coordinates {r,z,(f)}, we specify the injected poloidal 
flux function \E' as, (where \1' = rA^ and A^ is the component of vector potential), 

\l>(r, z) = r^ exp(-r^ - z^) , (10) 

up to a normalization coefficient Bq. Note that one can easily introduce another scale variable 
in the z— direction which separates the scaling between radial and vertical directions. Here, 
we have taken them to be same for simplicity. In the global distributions, we envision that 
the poloidal fiux and current will point along one direction in one region (say, for r < 1) and 
"return" in the opposite direction in another region (say, for r > 1), so the total net poloidal 
fiux and current are zero. 
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The above-mentioned key ingredients are motivated by the general framework of treating 
the black hole accretion disk system as a dynamo "machine", from which we use Eq.(lO) 
as a simple representation of the axisymmetric part of the poloidal flux that might come 
from the dynamo process. The continuous shearing of the poloidal flux lines by the disk's 
differential rotation generates the toroidal magnetic fields which correspond to a large scale 
poloidal current. Although the system scale we are modeling here is much larger than the 
actual size of the black hole accretion disk system, we envision that there still exists such 
global poloidal flux and current distributions on these scales. 

The poloidal fields, up to a normalization coefficient Bq, are: 

Anj,. = --^ = 2zr exp{-r^ - z^) , (11) 

r oz 

An,, = l^ = 2{l - r-) eM-r' - z') . (12) 

The poloidal fields have an 0-point at {r,z) = (1, 0) where both -Binj,r and -Binj,^ are zero. In 
specifying the toroidal component Anj.^, we adopt the approach of writing 

Anj,^ = fm/r , (13) 

where /(^E') is an arbitrary function of ^. Physically, this means that we have made the 
0— component of the Lorentz force (J x B)^ of the injected magnetic fields equal to zero. 
However, this does not prevent rotation being developed during the evolution by the Lorentz 
force of the combined fields and currents from both the injected fields and the already existing 
fields. Specifically, we choose /(^E') = tt^P, which means 

-Sinj,</> = tt^E'/'T = arexp(— r^ — z'^) , (14) 

where a is a constant and it has the units of inverse length scale. Other forms of /(^) are 
certainly possible, which will give different pitch profiles for the magnetic fields. 

Fig. 1. — left: Shown is the radial component of the J x B force as a function of radius r in 
the equatorial plane with a = 1, 3, and 5 for solid, dot-dashed and dashed lines, respectively. 
right: Shown is the vertical component of the J x B force as a function of radius r at the 
z = 0.5 plane with a = 1, 3, and 5 for solid, dot-dashed and dashed lines, respectively. 

This simple initial configuration has several properties. Firstly, the magnetic fields are 
"localized" because the fiux and energy drop exponentially as a function of distance from 
the center. This is an important difference from many previous studies where large scale 
background magnetic fields are often assumed to be present initially. In the poloidal plane, 
this global dipole-like configuration gives zero net fiux and current. 
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Secondly, the parameter a roughly specifies the fiux ratio between toroidal and poloidal 
components. The poloidal fiux between < r < 1 is 

$p = 27r^(l, 0) = 27r/e ^ 2.3 , (15) 

and the toroidal fiux threading through the {x, z} plane with a; > is 

$t = av^/2 ^ 0.9a . (16) 

So, the poloidal and toroidal fiuxes are roughly equal when a ~ 2.6. In astrophysical jets 
powered by AGNs, it is likely that magnetic fields are highly wound up by the disk rotation, 
this will imply that a >> 1. The poloidal current fiowing through the middle plane {z = 0) 
with < r < 1 is 

J^ = 27r/(^(1,0)) = 27ra/e^2.3a . (17) 

The total toroidal current is 

r r 

Jsdrdz = 2\fT\ . {li 




The ratio of poloidal (J^) to toroidal current (/^) is ay/n/e ~ 0.65q;, i.e., the poloidal current 
becomes dominant when a ^ 1. 

Thirdly, it is important to realize that this initial field is not in force equilibrium because 
J X B 7^ 0. The current density J = V x Bjnj is given as 

Jr = 2azr exp(— r^ — z'^) , (19) 

J, = 2a(l - r^) exp(-r2 - z'^) , (20) 

J0 = 2r(5 - 2z^ - 2r^) exp{-r^ - z^) . (21) 

Along the equatorial plane (2; = 0), the radial force density is given as 

F^(^ = 0) = 2r(l-r2)(10-a^-4r2)exp(-2r2) . (22) 

It can be seen that at large r, F^ > 0, meaning that the magnetic fields always expand 
outward due to the "hoop" forces. At small r and small a, the force is radially outward as 
well, whereas for small r but large a, the force is radially inward (i.e., "pinching"). Figure 
1 (left) illustrates these effects along the equatorial plane for different a values. Along the 
vertical direction, 

F, = 2zr^{a^-10 + 4:z'^ + 4:r^)exp{-2r^-2z'^) . (23) 



For a^ > 10, Fz is always positive (negative) for z > 0{z < 0), driving the fields away from 
the mid-plane. For smaller a, fields at large spherical radii will still expand away from the 
origin. Figure 1 (right) shows the vertical force as a function of r at 2; = 0.5 for different a. 

For the sake of completeness, we also give the total magnetic energy of this field: 

{5 + a^)n 



Eo = j{By2)dV = yi 



0.5(5 + a^) , (24) 



out of which the first and second terms are from the poloidal and toroidal magnetic field 
components, respectively. 

In summary, we have devised a magnetic field structure that is governed by the poloidal 
flux ^ and poloidal current J^. This magnetic field is localized in space, containing both 
poloidal and toroidal fluxes which are interwoven (i.e., finite helicity), and will dynamically 
evolve (i.e., not in force-free equilibrium), especially along the vertical direction in the large 
a limit. As we will show later in detail, it is this vertical expansion along the symmetry axis 
that forms the basis of modeling the large scale structure of jets. 



2.3. Units for 3D Simulations 

The initial magnetic field structure is embedded in a plasma background with finite gas 
pressure and density, both of which are taken to be uniform for simplicity in this study. We 
take Pg = Pq, p = po, and the sound speed Cg = \/^Pj p ~ 1.3a/-Po/Po- The magnetic 
field strength normalization parameter is B^. Simulations are performed in the Cartesian 
geometry with a uniform mesh. The simulation domain is a cube and the resolution is usually 
taken as 240 x 240 x 240 or higher. 

For the following runs, we typically choose po = 1, -Po = 1, and B^ = \. To put these 
numbers in an astrophysics context: suppose that the code unit po = 1 corresponds to a 
background density of 3 x 10^^ particles/cm^, and the background temperature is 7 keV. 
This means that the sound speed in the code unit Cs = ^/7 ~ 1.3 corresponds to 1.1 x 10^ 
cm/s. The combined density and temperature relation also means that the magnetic field 
-Bq = 1 in the code units corresponds to a physical magnetic field of ~ 20pG. Suppose the 
code unit L = 1 corresponds to 15 kpc, then the simulation duration t = 5 corresponds to 
~ 9.2 X 10^ yrs, and a total energy of 200 corresponds to ~ 6.2 x 10^^ ergs. The poloidal 
flux as given in Eq.(15) is ~ 9.4 x 10'^'^ Gcm^, and the poloidal current as given in Eq.(17) 
is ~ 1.7a X 10^^ Ampere. 
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Results 



We present the simulation results for two cases: the first has impulsively "injected" 
magnetic fields (i.e., as initial conditions) for a = 20. The second case is with a continuous 
injection of magnetic fields. 

Fig. 2. — The impulsive injection case with a = 20. (left): Shown is the azimuthally averaged 
pressure in the {r, z} plane. Azimuthally averaged poloidal flux surfaces (white-color) are 
overlaid as well. They are at t = 0,0.5, 1, respectively, (right): Same as in the left panel, 
except that t = 3,4,5, respectively. Note the cartesian computational box is [—12, 12] in all 
three directions for the simulation, and axis limits and colorbars have been changed between 
left and right panels. 



3.1. Impulsive Injection with an initial a = 20 

We choose Pq = po = Bq = 1 and a = 20 for this run and the initial B is of the form of 
Binj. The injection coefficients 7b and 7p are taken to be zero (i.e., no continuous injection 
of B or p). The computational box is taken to be [—12, 12] in all three directions. These 
parameters give the maximum initial magnetic field \B\ ^ 8.6 at (r, z) ^ (0.7, 0). 

Fig. 3. — The impulsive injection case with a = 20. Shown is the overall energy evolution of 
different components. The magnetic, kinetic and changes in thermal energies are shown as 
solid, dashed, dot-dashed curves, respectively. The total energy is shown by the thin solid 
curve at the top. As the system evolves, the initial magnetic energy is converted to both the 
thermal energy and the kinetic energy of the plasma. 

The "free" energy of the system is initially all magnetic. Since the initial magnetic field is 
out of force equilibrium (see Figure 1), evolution is expected. Because the evolution remains 
quite axisymmetric (see below), in Figure 2, we show the time sequence of azimuthally 
averaged pressure distribution, overlaid with azimuthally averaged poloidal fiux function 
\E'. In Figure 3, we show the overall energy evolution for different components. The initial 
magnetic energy is given by Eq.(24) with a = 20. The kinetic energy is pv'^/2 integrated 
over the whole volume, and the change in thermal energy is [Pg{t) — -Pg(0)]/(7 — 1) integrated 
over the whole volume, where -Pg(t) and -Pg(O) are the pressure distributions at time t and 
t = 0. 

The whole evolution can be approximately separated into two stages. The first stage 
(0 < t < 1) is a rapid conversion of ~ 70% of the available magnetic energy through the 
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work done by the J x B force on the surrounding plasmas, which pinches all the flux lines 
within r = 1 radially inward and expands flux lines vertically away from the mid-plane. 
Since a in this case is quite high, the J x B force is consequently so strong that it initially 
drives both a radially converging shock and an axially expanding shock, as indicated by the 
central red region and the red region surrounding the magnetic fields in the left panel of 
Figure 2. These shocks heat the plasmas efficiently and generate large amounts of entropy. 
The second stage is much slower, represented by the slow dissipation of magnetic energy and 
further gradual heating of the plasma. The magnetic field expansion gradually slows down. 
In the meantime, a compressional shock wave, shown as the outer red region in the right 
panel of Figure 2, continues with its expansion. Figures 4 and 5 show the time evolution of 
the azimuthally averaged radial profiles of the radial velocity Vr and the plasma pressure Pg. 
A shock front, generated from the initial impulsive release of the magnetic energy, can indeed 
be seen. This shock heats the plasmas surrounding the magnetic structure as indicated by 
the propagating pressure front. The continued expansion, however, also makes a central, low 
pressure cavity. Note that the gas pressure changes from this wave propagation are quite 
strong (~20-40%). 

Fig. 4. — The impulsive injection case with a = 20. Shown is the radial distribution of the 
azimuthally averaged (cylindrical) radial velocity f ^ at t = 0.5, 1.5, ...,4.5, respectively. The 
profiles are taken at z = 2.4. The front generated by the sudden release of magnetic energy 
is moving supersonically. 

Fig. 5. — The impulsive injection case with a = 20. Shown is the radial distribution of the 
azimuthally averaged pressure at t = 0.5, 1.5, ...,4.5, respectively. The profiles are taken at 
z = 2.4. At early time, the "pinching" of inner magnetic fields greatly increases the pressure 
near the central axis. This is followed by expansion which creates a central low pressure 
(and density) region. The front is moving supersonically. 

Fig. 6. — The impulsive injection case with a = 20. Shown is the magnetic fields \B\, density 
and pressure distributions in the {x, z} plane (y = 0) at t = 5. The overall structure remains 
axisymmetric. The magnetic field expansion creates a central low density and low pressure 
region. 

Figure 6 shows the distribution of the magnitude of magnetic fields \B\, density p and 
gas pressure Pg at t = 5 in the {x,z} plane {y = 0). These distributions suggest that the 
final state is still quite axisymmetric. The magnetic fields form an elongated structure due 
to both the pinch at the central region and the vertical expansion along the symmetry axis. 
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Fig. 7. — The impulsive injection case with a = 20. Shown is one magnetic field line at 
t = 5. The field line goes up in a tightly wound central helix and comes back in a loosely 
wound helix. 

Fig. 8. — The impulsive injection case with a = 20. [Upper): The radial distribution of 
magnetic field components B^ and B^. (Lower): The radial distribution of the plasma 
pressure Pg, the magnetic pressure Pb, and the total pressure Ptot = Pg + Pb- Both panels 
are taken at z = 3 and t = 5. 

Figure 7 depicts a 3D field line, showing its going up in a central tightly wound helix and 
coming down outside in a loosely wound helix. Based on the rate of field twists, one might 
expect that this helix should be kink unstable. The hollow pressure profile as shown in 
Figure 6, however, may be a stabilizing infiuence. 

Figure 8 shows the radial distribution of azimuthally averaged magnetic field compo- 
nents (top) and different pressure components (bottom) at 2; = 3 and t = 5. The expansion 
of the magnetic fields from an initial sphere with a radius of ~ 1 to the elongated structure 
have created a very low density cavity since plasmas are tied to the expanding magnetic 
fields. The pressure within this cavity is mostly from the shock heating. The cavity is dom- 
inated by magnetic pressure (the plasma {3 = Pg/Ps is ~ 10% near the axis). A careful 
force balance analysis shows that the central region is nearly force free (more detailed results 
will be presented elsewhere). Going to large radii (r ~ 1 — 1.5), however, the structure is 
confined by the gas pressure, which prevents further expansion of the magnetic fields. 

Since the density is quite low inside the cavity, there is some numerical heating, which 
amounts to an error of ~ 2% increase in the total energy (the top thin solid curve in Figure 
3). Some of the poloidal fiux lines at late times (the plot at t = 5 of Figure 2) indicate 
numerical magnetic reconnection at the mid-plane. 



3.2. Continuous Injection 

Fig. 9. — The continuous injection case with an injection a = 15. Shown is the evolution 
of azimuthally averaged poloidal fiux lines during and after the magnetic field injection. 
Different plots are for t = 0,2.5,4,5.5, and 7. The injection is terminated at t = 2.5. Each 
subplot has ten, evenly spaced contour lines. 

For this case, we adopt the functional form of Bjnj for both the initial and injected 
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magnetic fields. For the initial fields, we choose Pq = po = 1, Bq = 1, and a = 3. For the 
continuous field injection, we use Bq = 0.2, a = 15, 7b = 3, 7p = 0.1, r^ = 0.5, and t^^j = 2.5. 
The computational box is [—12, 12] in all three dimensions. This run uses a grid of 320^. To 
put these numbers in an astrophysics context, using the conversions given in the beginning 
of this section, we have: the simulation's tinj = 2.5 corresponds to ~ 4.6 x 10'^ yrs, and a 
total energy of 200 corresponds to ~ 6.2 x 10^^ ergs. The mass injection rate 0.1 corresponds 
to ~ IMq/ut and the magnetic field Bq = 0.2 corresponds to ~ 4/xG. All these parameters 
are roughly consistent with the physical conditions in the central region of galaxy clusters. 
However, note that the background medium is not stratified in the current simulation, so the 
total thermal energy within a computational volume of 10^ would correspond to ~ 3 x 10^° 
ergs. In other words, the magnetic fields are expanding into a volume with a high thermal 
pressure background. 

Figure 9 shows the evolution of azimuthally averaged poloidal flux lines in the {r, z} 
plane. The initial poloidal flux has a maximum of ~ 0.37, then the poloidal flux increases to ~ 
0.9 at t = 2.5 when the injection finishes. Magnetic fields continue to expand somewhat after 
the injection has stopped. Figure 10 shows the energy evolution for different components. 
The injected energy is essentially all in magnetic fields, with a very small amount in the form 
of thermal and kinetic energies of the injected plasmas associated with p^^y There is again 
a little numerical heating at the late time (see the top two curves in Figure 10) due to the 
existence of extremely low density regions. 

Part of the injected magnetic energy is gradually converted into both the kinetic energy 
and the changes in the thermal energy of the surrounding plasmas. The system is approaching 
a quasi steady state with very slow or no changes in different energy components. We 
interpret this as that the injected magnetic fields run out of their "drive" since the magnetic 
field injection stops at t = 2.5. Subsequently, the magnetic fields have come to be in rough 
force balance with the surroundings. Of course the whole system can not be steady, it is still 
expanding, though very slowly. Figure 11 shows the distribution of magnetic fields, density 
and gas pressure at t = 7. A magnetically dominated cavity is evident. 

Figure 12 shows the radial distribution of azimuthally averaged magnetic field com- 
ponents and pressures at t = 7 and z = 3. Again, we have a central column which is 

Fig. 10. — The continuous injection case with an injection a = 15. Shown is the evolution of 
different energy components. The thicker solid, dashed and dot-dashed lines are for magnetic, 
kinetic and changes in thermal energies, respectively. The thinner dashed and solid curves 
at the top are for the total injected energy and the total energy in the computational box, 
respectively. 
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magnetically dominated, with nearly the force-free condition. At a larger radius (r ~ 1.5), 
a small gas pressure gradient is present to balance the residual magnetic pressure gradient. 
This hollow pressure profile may be helpful to stabilize the whole structure. 

It is interesting to note that more than half of the injected magnetic energy still remains 
in the system for this set of parameters. Compared with the impulsive injection case with 
a = 20, the two systems have about the same initial "free" energy, and we can see that the 
continuous injection case retains a much larger fraction of the magnetic energy and gives a 
slightly larger magnetically dominated volume. 



4. Discussions 

In this paper, we have presented a method of injecting free energy in the purely magnetic 
form. The non-force-free nature of these magnetic fields causes expansion in a self-pinched, 
collimated fashion. In relating to astrophysical jets, we assume that the whole radio galaxy 
can be divided into: a central region where the black hole accretion system resides, a prop- 
agation region which contains the jets and outflows supplied by the central region, and a 
termination region where the lobes are formed. Based on this assumption, we suggest that 
the approach taken here could be used to study both the propagation and termination re- 
gions. We have especially concentrated on the magnetically dominated regime, different from 
the previous studies where energy is predominantly carried by the kinetic energy of the flow. 

The 3D MHD simulation results for a few relatively simple cases are presented here, 
mostly for the purpose of illustrating the salient features of the particular functional form of 
the injected magnetic flelds and how the subsequent evolution might be related to it. As such, 
we have only explored a limited set of parameters and their influence on the magnetic fleld 
evolution. Many important issues are not addressed here but will be presented in forthcoming 
papers. For example, it will be important to know the stability of such magnetic structures. 
The fleld line plots show that they are highly wounded, suggesting kink unstable, yet the flnal 
distributions seem to suggest that they remain axisymmetric and stable. The hollow pressure 
proflle might have a stabilizing influence, but more detailed studies are required in order to 

Fig. 11. — The continuous injection case with an injection a = 15. Shown is the magnetic 
field, density and pressure distributions (from left to right) at t = 7. The central low density 
and pressure volume is occupied with and dominated by magnetic fields. A compressional 
shock wave generated during the magnetic field injection stage has propagated away from 
the central, magnetized region. 



-14- 



draw firm conclusions on the stability. In addition, how background environment can affect 
the propagation speed and general morphologies, how the final configuration depends on the 
injection rates, and what determines the total amount of magnetic energy dissipation, etc. 
It is thus premature to compare the current simulation results directly with observations of 
radio jets/lobes, though this is our ultimate goal. 

HL acknowledges useful discussions with D. Ryu and M. Nakamura. This research 
was performed under the auspices of the Department of Energy. It was supported by the 
Laboratory Directed Research and Development Program at LANL and by IGPP at LANL. 
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